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ABSTRACT 


^  A  computer  code  has  been  developed  for  solving 
incompressible  two-dimensional  axisymmetric  time-dependent 
viscous  fluid  flow  problems  involving  up  to  two  free 
surfaces.  Heuristic  models  for  turbulence  are  employed 
to  extend  the  method  to  indefinitely  high  Reynolds 
number.  Scalar  quantities  (heat  and  solute  concentrations) 
are  also  followed,  and  the  fluid  may  be  slightly  non- 
homogeneous  in  the  Boussinesq  approximation.  The  method 
is  a  second-order  space,  forward  time  explicit  finite- 
difference  scheme.  Free  surfaces  are  treated  using  the 
MAC  ("Marker-and-Cell")  technique. 
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1 .  INTRODUCTION 

The  MACYL6  computer  program  was  created  in  response  to  a 
requirement  to  study  the  events  following  the  detonation  of 
a  nuclear  device  far  below  the  surface  of  the  sea.  Of  principal 
interest  are  the  pulsation  and  upward  migration  under  gravity 
of  the  steam  bubble  generated  by  the  explosion  and  the  resulting 
redistribution  and  upward  translation  of  the  radioactive  bomb 
debris.  The  MACYL6  code,  while  adequate  to  treat  the  explosion 
problem,  is  also  applicable  to  a  wide  variety  of  applied 
problems  in  various  areas.  It  is  not  the  purpose  of  the 
present  paper  to  discuss  such  applications,  but  rather  to 
describe  the  numerical  method  itself  in  detail.  Applications 
to  particular  cases  (including,  of  course,  the  explosion  problem) 
will  be  published  in  subsequent  reports. 

The  MACYL6  code  was  developed  in  an  evolutionary  fashion; 
an  earlier  and  more  primitive  version  is  described  in  Pritchett 
(1967).  In  earlier  forms,  it  has  proved  useful  in  the  calcula¬ 
tion  of  explosion  phenomenology  and  has  Leen  successfully 
compared  with  measured  field  data  (Pritchett  and  Pestaner,  1969). 
The  latest  version  is  a  great  deal  more  complicated  than 
those  published  earlier,  differing  from  them  in  two  essential 
respects.  First,  scalar  transport  equations  are  introduced 
(in  addition  to  the  equations  of  fluid  motion)  which  permit 
the  simultaneous  calculation  of  the  space-time  distributions 
of  such  quantities  as  temperature,  salinity,  bomb-debris 
concentration,  and  the  like.  The  effects  of  these  distribu¬ 
tions  may  then  be  reintroduced  into  the  fluid  mechanics  by 
allowing  the  fluid  density  to  be  slightly  dependent  on,  say, 
temperature  and  salinity.  This  permits  such  effects  as 
oceanic  stability  and  thermohaline  convection  to  be  included 
in  the  scheme.  Second,  and  more  important,  a  recently-developed 
heuristic  model  for  fluid  turbulence  is  an  essential  part  of  the 
method.  This  model  simulates  the  effects  of  turbulence  on 
the  flow  through  distributions  of  turbulent  energy  and 


scales  of  turbulence.  The  essential  features  of  the 
model  are  summarized  in  section  II;  for  justification, 
details,  and  comparisons  with  measured  data  the  reader 
is  referred  to  the  original  papers  (Gawain  and  Pritchett, 
1969,  1970). 

The  numerical  method  is  an  explicit  forward-time 
finite-difference  scheme.  The  momentum  equation  uses 
a  nine-point  second-order  space  difference  representa¬ 
tion  of  the  advection  terms,  and  the  general  scalar 
transport  equation  uses  the  more  stable,  but  somewhat 
less  precise  "upstream"  or  "donor-cell"  method.  These 
are  described  at  length  in  section  III.  The  pressure- 
velocity  (rather  than  the  stream  function-vorticity ) 
formulation  of  the  momentum  equation  is  used.  The  MAC 
( "Marker-and-Cell" )  free-surface  treatment  is  employed; 
this  scheme,  developed  at  Los  Alamos  (Welch,  Harlow, 
Shannon  and  Daly,  1966) ,  represents  tha  fluid  by  a 
number  of  massless  "marker  particles"  which  move  with 
the  flow  through  the  Eulerian  mesh  and  thereby  specify 
the  position  of  the  free  surface.  The  treatment  is  in 
two-dimensional  axisymmetric  cylindrical  coordinates, 
and  the  flow  is  assumed  to  be  truly  two-dimensional 
(that  is,  without  swirl) .  The  computer  program  itself 
is  largely  written  in  FORTRAN  IV,  with  some  portions  in 
assembler  language  to  increase  speed;  it  is  currently 
operating  on  Control  Data  6600  equipment. 


II. 


THE  GOVERNING  EQUATIONS 


a.  The  Navier-Stokes  Equations 

In  Eulerian  coordinates  (that  is,  coordinates  fixed  in 
space  rather  than  moving  with  the  flow)  the  equations 
of  motion  for  an  incompressible,  homogeneous  viscous  fluid 
in  a  gravitational  field  may  be  written  in  vector  notation 
as  follows: 

V  •  u'  =  0  (II-l) 

+  v‘  (U'U')  =  -V4>  '  +  uV2u'  +  g  (II-2) 

where  U '  =  velocity 

$ '  -  pressure/density  (density  being  constant) 
u  =  fluid  kinematic  viscosity 
g  =  acceleration  of  gravity 

These  equations  express  the  principles  of  mass  and  momentum 
conservation,  respectively.  The  two  terms  on  the  left  of 
(II-2)  represent  the  total  time  rate  of  change  of  momentum 
for  an  element  of  fluid  moving  with  the  flow.  The  first 
term  on  the  right  is  the  rate  of  momentum  production  due  to 
normal  pressure  forces,  the  second  expresses  the  rate  of  momentum 
diffusion  by  viscosity,  and  the  third  is  the  rate  of  momentum 
production  by  gravitational  forces. 

b.  Scalar  Transport  Equations 

Other  transport  equations  which  will  prove  useful  later 
in  the  development  are  the  transport  equations  for  heat  and 
solutes : 
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where  T  '  •»  temperature 

S  '  -  solute  concentration 
ic  -  thermal  dlffusivlty 
C  ■  molecular  diffusion  coefficient 
nT  ■  rate  of  temperature  change  due  to  heat  sources 
The  last  term  in  equation  (11-3)  (the  source  term)  will 
not  be  pursued  further  at  this  time;  it  represents  the 
rate  of  dissipation  of  kinetic  energy  to  heat. 

c.  Turbulence  and  the  Reynolds  Stresses 

In  incompressible  flow,  the  equations  of  continuity 
and  momentum  (II-l  and  II-2) ,  along  with  the  boundary 
conditions,  in  principle  establish  completely  the  entire 
fluid  motion.  If  the  flow  is  turbulent,  however,  the 
detailed  motion,  although  theoretically  determinate, 
becomes  so  complex  that  its  actual  calculation  would 
involve  an  overwhelming  amount  of  computation.  Further¬ 
more,  in  general  the  results  of  interest  are  certain  average 
properties  of  the  flow,  and  the  large  mass  of  additional 
detailed  information  available  is  usually  neither  required 
nor  desired. 

In  Cartesian  tensor  notation,  equations  (II-l)  and 
(II-2)  may  be  written: 


2-  (u  r  )  =  o 

8x.  1 
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(II-5) 
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(II-6) 


If  we  express  the  total  velocity  and  pressure  fields  as 
the  sum  of  their  mean  and  fluctuating  components, 


u.'  =  u. 

1  1 


+  u. 


(II-7 ) 


I*  *  #  *  t 


(Total  •  ooaii  «  turbulent) 

substitute  these  into  (II-))  and  (It-))  ant  ensemble- 
average  the  results*  we  obtain  continuity  and  nowentun 
equation*  for  the  nean  flout 

J  (Ut)  *  0  (Il-tl 


JUt  a  J  JOj  JU* 

TT“  ♦  lx”  <uiV  *  "  llj  *  V 

(11-10) 

(overbara  denote  ensemble-averaged  quant  it  let) 


Further  manipulation  yields  an  energy  equation  for  the 
turbulent  f luctuationsi 


It  will  be  noticed  that  the  Man  flow  pone n tun  equation 
(11-10)  resembles  the  original  total  none n tun  equation  (II-)) 
except  for  the  presence  of  the  unknown  Reynold*  stresses 
(-opr  ) .  it  is  customary  and  appropriate  to  postulate  that 
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these  Reynolds  stresses  can  be  adequately  related  to  the 
mean  flow  strain  rates  through  the  law: 


'uiuj  ■  "  j  Vk  sij  +c(^j +  (I 

where  ■  0  for  i  ^  j  and  ■  1  for  i  «  j;  e  is  the 

so-called  eddy  kinematic  viscosity.  The  mean  flow 
momentum  equation  is  then: 


sui  ♦ 
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(11-13) 


The  above  averaging  process,  as  has  been  seen,  results 
in  a  great  simplification  of  the  physical  problem  ,  but 
also  involves  a  significant  and  irretrievable  lors  of 
essential  information.  This  is  apparent  in  the  appearance 
of  the  unknown  Reynolds  stresses,  or  alternatively,  of 
the  unknown  eddy  viscosity  distribution.  To  define  determi¬ 
nate  solutions,  additional  relations  are  needed.  Regrettably, 
such  supplementary  relations  cannot  be  established  from 
the  original  equations  (II-l)  and  (II-2)  by  any  purely 
deductive  process,  and  therefore  empirical  hypotheses  are 
an  unavoidable  necessity.  From  another  point  of  view, 
the  averaged  equations  of  motion  show  the  effect  of  the 
Reynolds  stresses  on  the  mean  flow,  but  the  reciprocal 
effect  of  the  mean  flow  on  the  Reynolds  stresses  is  lost 
in  the  averaging  process.  Hence,  some  adequate  hypothesis 
must  be  found  to  approximate  this  relationship. 

For  this  purpose,  a  recently  developed  model  for  turbu¬ 
lent  flow  (Gawain  and  Pritchett,  1969;  1970)  is  used  in 


the  present  work.  It  is  similar  in  concept  to  earlier 
work  of  Prandtl  (1945)  and  is  somewhat  similar  to  a  scheme 
proposed  at  Los  Alamos  (Harlow  and  Nakayama,  1967;  Hirt ,  1968). 
Only  the  barest  essentials  of  the  scheme  will  be  sketched 
here;  for  physical  rationale  and  other  details  the  reader 
is  referred  to  the  original  papers. 

The  basic  postulate  in  the  model  is  that  the  eddy 
viscosity  assumption  (11-12)  is  appropriate ,  and  that  the 
eddy  viscosity  may  be  adequately  represented  by  a  relation 
of  the  form: 

e  =  «  /{Tin  A  =  af~TE~  A  (11-15) 

where  a  is  a  slowly  varying  dimensionless  function  which 
may  be  taken  as  a  constant  outside  boundary  layer  regions. 

A  is  a  "scale  of  size"  associated  with  the  mean  flow  which 
may  vary  from  point  to  point  in  space  and  time  and  which 
will  be  discussed  later. 

Combination  of  the  basic  turbulent  energy  equation 
(11-11)  with  the  eddy  viscosity  postulate  (11-12)  yields: 


(11-16) 
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(mean  flow  strain  rate) 
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where 


(11-18) 
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(mean  turbulent  kinetic  energy)  (11-19) 


The  terms  on  the  right  of  this  energy  equation  represent, 
respectively,  turbulent  energy  production  corresponding  to 
the  work  done  by  the  mean  flow  against  the  Reynolds  stresses, 
dissipation  of  turbulent  energy  to  heat,  "turbulent  diffusion" 
of  energy,  and  molecular  diffusion.  In  the  heuristic  model, 
the  last  term  is  neglected,  as  it  is  vanishingly  small  at 
high  Reynolds  numbers.  The  remaining  terms  are  approximated 
as  follows: 
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(production) 


-  u  —  (dissipation) 

X2 

+  (ye  '!“•]  (diffusion) 

(11-20) 


where  X  is  the  so-called  "dissipation  length"  or  "turbulent 
microscale",  and  y  is  another  slowly  varying  function  of 
the  same  sort  as  a  (see  equation  11-15)  which  is  constant 
outside  boundary  layers. 

To  complete  the  model,  it  is  now  only  necessary  to 
establish  the  "macroscale",  A  (equation  11-15)  and  the 
"microscale",  X  .  We  define  a  generalized  mean  flow  strain 
rate  ft  and  a  generalized  mean  flow  strain  rate  gradient  ft  '  as 
follows: 


ft2  = 


1 

I 


ij 


(see  equation  11-18) 


(II-21) 


-8- 


(11-22) 
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From  these  definitions ,  the  following  quantity  can  be  obtained: 
m*  -  i  iioiwio* 


W> 


T  <TT>  (y^ 

xi  i 


We  now  define  the  macroscale  as  follows: 


A2 (x,t) 


I2 (x,t) 
J  (x,t) 


where 


1 2  (x ,  t )  =  J  w(x,x*  ,t)  R4  (x'  ,t)  dv' 
all  space 

J2(x,t)  =  j  w(x,x',t)  (RR'  (x' ,t) ) 2dv' 
all  space 


and 


w(x,x',t)  = 


exp 


_  f(x-x')Mx-x')  ) 

1  A2 (x,  t)  J 


exp  ]dv. 

t  A2  (x» t) 


all  space 


(11-23) 


(11-24) 


(11-25) 


(11-26) 


(11-27) 


Difficulty  will  of  course  be  experienced  in  calculating  the 
distribution  directly  from  an  arbitrary  velocity  distribution, 
since  A  appears  on  both  sides  of  the  defining  equation.  To 
calculate  A  explicitly,  the  following  iteration  process  is  used. 
Let  the  (n+l)-th  approximation  to  A  be  defined  as: 


A 


2 

n+1 


(x,t) 


(x,t) 

(x,t) 


(11-28) 


where 


In+1  (*»t>  "  j  wn(x,x' ,t) n4(x' ,t)  dv' 
all  space 

Jn+l(*'t}  =  j  wn(x,x',t)<nn'  (x',t))2dv' 
all  space 


(11-29) 


(11-30) 


and 


exp 


wn(x,x' ,t)  = 


(x-x* ) • (X-X* )  ) 
A* (x,t) 


exp 


L  (x-xMlK-x1)  1 


all  space 


A*  (x,t) 


dv1 


(11-31) 


To  start  the  iteration  process,  we  take  Ao®"0  everywhere. 
Then  in  computing  Ai,  we  find  that  the  weighting  function 
w0  is  the  same  everywhere,  so  we  obtain  simply: 

j  ft4dv 

^  all  space _ 

(nn')2dv  di-32) 

all  space 

a  constant  independent  of  position.  Thus,A2  is  the  first 
non-constant  approximation  to  A  that  we  obtain.  As  was 
discussed  in  the  original  paper,  the  rapid  convergence  of 
the  An's  justifies  the  approximation 


A (x, t) *  A2 (x, t) 


(11-33) 


This  procedure  is  also  used  in  the  present  work. 
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The  dissipation  length  X  is  formulated  as  follows. 
Consider  two  lengths  and  L2  defined  as: 


2eT 


2E 

J 


L2  = 

where  J  is  as  defined  above. 


(11-34) 


(11-35) 
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We  now  form  the  relation: 


(11-36) 


where  8  is  yet  another  slowly  varying  dimensionless 
function  which  reverts  to  a  constant  outside  boundary 
layers.  Numerical  values  that  best  match  experimental 
data  for  regions  outside  the  turbulent  boundary  layer  are: 

(11-37) 


a  =  0.065 
3.7 
1.4 


6 

Y 


(11-38) 

(11-39) 


d.  Turbulent  Scalar  Transport 

As  has  been  seen,  the  scalar  transport  equations  for 
heat  and  solutes  are  of  the  form: 


(U'.Q') 


c  92Q' 

SxTSxJ 


+  n 


Q' 


(11-40) 


where  Q'  is  the  instantaneous  scalar  concentration,  the  U.' 
are  the  instantaneous  velocity  components,  C  is  a  molecular 
diffusivity,  and  n^,  is  a  source  or  sink  function.  These 
may  also  be  expressed  as  the  sum  of  mean  and  fluctuating 
components : 


0'  = 

Q  +  q 

(11-41) 
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As  was  done  for  the  momentum  equation,  the  scalar  transport 
equation  may  be  averaged.  The  result  is: 


(UjQ) 


+  IL 


In  analogy  to  (11-12)  we  postulate: 


(11-43) 


(11-44) 


which  defines  the  turbulent  diffusion  coefficient  D^. 

Although  the  molecular  dif fusivities  of  the  various 
scalar  quantities  (such  as  heat,  salt,  etc.)  may  be  quite 
different,  the  turbulent  diffusion  process  is  a  property 
of  the  turbulence  field,  rather  than  of  the  particular 
substance  being  diffused.  Furthermore,  at  high  Reynolds 
number,  it  is  to  be  expected  that  the  turbulent  diffusion 
coefficient  will  be  vastly  larger  than  the  molecular 
diffusion  coefficient.  Thus,  for  our  purpose,  it  is 
sufficient  to  write: 

(u.q)  =  ,JL.  (d  )  +  n 

3t  3x7  vuju'  3x7  'UT  3x7'  Q  (11-45) 

We  have  already  approximated  the  turbulent  diffusion 
coefficient  for  turbulent  kinetic  energy  as  ye,  where 
Y  =  1.4.  Numerous  measurements  (Corrsin  and  Uberoi ,  1947; 
Schlicting,  1960)  ,  suggest  that  the  "turbulent  Prandtl  number" 
(turbulent  thermal  dif fusivity/eddy  viscosity)  generally 
takes  on  values  between  1  and  2  for  most  flow  situations. 

This  suggests  that  it  is  appropriate  to  assume  that 


Dt  =  ye  =  1.4e 


(11-46) 
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that  is,  that  the  turbulent  diffusivity  for  turbulent  energy 
applies  to  all  scalar  turbulent  diffusion.  Thus  finally  we 
obtain: 


II  +  HT  (ujQ>  '  +  "q  (I1-47 

The  source  term  for  turbulent  energy  has  already  been 
defined  -  it  is  simply  the  work  input  from  the  mean  flow  minus 
the  rate  of  dissipation  to  heat.  The  source  term  for  the  heat 
equation  reflects  the  dissipation  of  both  turbulent  energy 
and  mean  flow  kinetic  energy  (the  latter  is  usually  small) . 
Thus,  for  heat  transport,  we  write: 


3T  _3_ 

at  axj 


<ujT»  -  dr  <*4c7> 


♦  I  <«! 


2E  . 
X2-  } 


where  T  =  temperature 

o  =  specific  heat  of  fluid 


(11-48) 


e.  The  Boussinesq  Approximation 

So  far,  we  have  discussed  an  homogeneous  fluid,  that  is, 
a  fluid  in  which  the  molecular  viscosity  and  the  mass  density 
are  constants.  In  oceanic  problems  however,  for  example, 
even  a  slight  non-homogeneity  in  density  may  have  profound 
effects  upon  the  resulting  flow.  Although  the  density 
difference  is  quite  small,  the  interaction  of  the  gravitational 
acceleration  with  a  slight  density  stratification  may  set  up 
regions  of  stability  or  instability  which  are  important  in 
the  problem.  It  should  be  pointed  out,  however,  that  density 
fluctuations  caused  by  pressure  fluctuations  (i.e.  slightly 
compressible  behavior)  are  irrelevant  in  this  regard.  Only 
density  variations  due  to  intrinsic  properties  of  the  fluid 
parcel  (such  as  temperature  or  salinity)  have  such  a  stabilizing 
or  de-stabilizing  effect.  Therefore,  rather  than  actual 
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"in-situ"  density,  we  will  be  concerned  with  "potential 
density":  that  is,  density  referred  to  some  standard  pressure. 

Thus,  we  may  re-write  the  mean-flow  momentum  equation 
in  the  Boussinesq  approximation  (retaining  the  density  effect 
in  the  gravity  term  but  not  in  the  inertial  terms) : 


3U. 
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3 


(W 


3P 

**I 
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3Xj 


t (u+e) 


3U . 

‘*5 


3U. 

J^)] 


+  Kg* 


(11-49) 

where 

C  -  1  +  ^  -  f(T,S),  S  =  salinity  (11-50) 

P  0 

The  salinity  transport  equation  is  simply  (see  11-47) 

H  +  §57  <V>*  ffj>  <»-“> 

For  seawater,  for  example,  the  dependence  of  potential 
density  on  salinity  and  temperature  is  well  known  experimentally 
(Knudson,  1901) .  For  temperature  in  centigrade  degrees  and 
salinity  in  parts  per  thousand,  over  a  normal  range  of  conditions 
at  one  atmosphere  the  data  is  well  represented  by: 


*  innn  -  (T-3.98) 2  (T+283) 
p (kg/m  )  1000  503.57  (T+67.26) 


.S-0.03, 


,S— 0 . 03t 2 


+  (0.0634  +  1.4708  (f-^---^)  -  0 . 00157  (j-jfcjp) 

+  0. 0000398 (-fTyof) *  I x (— -r(18. 03-0 . 8164T+0 . 01667T2 ) 
x  (-0.2014+1.4708  (S~^°3)-0.00157  (S~°gj}3)  2 
+  0.0000398  (^-°^-g)  »)  +  .!  -  j^qq'(4 . 7867-0 . 098185T 


+  0 . 0010843T2 ) 3 


(11-52) 
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f .  Summary  of  Principal  Equations 


In  two-dimensional  axisymmetric  cylindrical  coordinates, 
the  field  equations  may  be  written  as  follows: 

Continuity: 


?  3F  (ru)  +  17 

Radial  momentum: 


8v 


=  0 


(11-53) 


)u  .  1  3  ,  2».  9,  » 

St  +  F  IF  (ru  >  +  I?  (uv)  * 


's  r  1  9  /  ^  9  U  v  s  ,  9  ( r-  it  ^  \  1 

21  F  IF  (rE  IF1"  ~r  *  IF  U  IF  1 


+  [£*(^ 
3z  1  '  3z 


3P 

3r 


3v.  . 

17)  1 


(11-54) 


Vertical  momentum: 


If  ♦  F  IF  (ruv)  +  IF'*2’  - 


2  I?  If  <"*  Si»  +  it<HiF>> 


+  r  IF  <§"  5F>1 


9P  4.  f„ 

IF  +  £s 


(11-55) 


Turbulent  energy: 


(ruE)  +  —  (vE)  =  -  (rve~)  +  —  (ye—  ) 
3 1  r  jr  gz  '  '  r  gr  '  TEgr'  gz  vyE3z  ' 


_  r,  2  2E 

+  efi2-  u— 
X2 


(11-56) 
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Heat: 


3T  ,  1  3 

7t  +  r  3r 


(ruT) 


1  3_ 
r  3r 


(rye 


3T 

r1 


/  3T. 


£  (21  +n! 
0  X2 


Salinity  and/or  other  solutes: 


(11-57) 


7E  +  r  Ir  (ruS)  +  h  (vS) 


13  t  3S\  3  t  3  S  * 

r  Tr  (rYe  Sr*  +  7?  (y£3F* 


(11-58) 


where : 

r  =  radial  coordinate 
z  =  vertical  coordinate 
u  =  radial  mean  flow  velocity 
v  =  vertical  mean  flow  velocity 
e*  =  e  +  u 

o  *  molecular  kinematic  viscosity 
P  =  total  pressure 

E  =  turbulent  kinetic  energy  per  unit  mass 
T  *  temperature 

S  *  salinity  (or  other  solute  concentration) 

£  =  density  deviation  parameter  =  (1  +  ^  )  =  prescribed 
function  of  S  and  T.  Po 

o  *  specific  heat 

e  =  eddy  viscosity  *  a. /Tit  A 

A  *  turbulent  macroscale  =  I2/J2 

ft2  =  (generalized  mean  flow  shear  rate)2  = 


•  tt'UMiH 


1 

m')2  -  \  uf^)1  +  <|§"> * > 

(11-61) 

I 

I2  (r ,  *)  -  //w(r>r>ti>i>)(n>(r,>r)),dr*di> 

(11-62) 

1 

//w(r,r' ,*#*')  dr'dx* 

I 

TP- 

aMr.z,  - 

(11-63) 

I 

^2  _  //n1*  r  dr  dx 

[ 

0  //  (nn') Jr  dr  dx 

(11-64) 

r 

X2  »  (turbulent  microscale)1  ■  u  -y 

6 (2EJ* )' * 

(11-65) 

1. 

[ 

(11-66) 

,  /rv7"  -f <*-*')*  ♦  <*-*')? 

w(r,r'(z,z')  -  f  (^-)  e  1  A»  1 

[ 

»(x)  -  25  f  e  -J*’  ‘l-'O*  •»  40 

/v  * 

0 

(11-67) 

I 

a  *=  0.065 

(11-68) 

I 

1  ,  , 

6  -  3.7 

(11-69) 

III.  THE  FINITE-DIFFERENCE  METHOD 


a.  The  Computing  Mesh 

Consider  the  region  of  interest  to  be  divided  into  a 

number  of  toroidal  cells  of  rectangular  cross-section  and 

not-necessarily-equal  cross-sectional  area  (see  Fig.  1) . 

If  we  denote  a  cell  by  the  indices  (i,j)  where  i  varies  with 

radius  and  j  varies  with  height,  we  may  establish  the 

following  nomenclature: 

Ar.  =  radial  dimension  of  cell  ij 

Az  .  =  vertical  dimension  of  cell  ij 

r^  =  distance  from  axis  to  center  of  cell  ij 

z  .  =  distance  from  bottom  of  mesh  to  center  of  cell  i j 
3 

r^_jy2  =  distance  from  axis  to  inner  boundary  of  cell  ij 

ri+l/2  =  distance  from  axis  to  outer  boundary  of  cell  ij 

z._jy2  =  distance  from  bottom  of  mesh  to  lower  boundary  of  cell  ij 

2 j+1/2  =  di-stance  froin  bottom  of  mesh  to  upper  boundary  of  cell  ij 

Ar^_jy2  =  distance  between  centers  of  cell  ij  and  cell  i-1  j 

Ar^+]y2  =  distance  between  centers  of  cell  ij  and  cell  i+1  j 

Az  .  ,  =  distance  between  centers  of  cell  ij  and  cell  i  j-1 

AZj+]y2  =  distance  between  centers  of  cell  ij  and  cell  i  j+1 

b.  Continuity  and  Momentum 

The  field  variables  are  calculated  at  specified  points 
within  the  cell  matrix  as  shown  in  Figure  2.  Scalar  quantities 
are  defined  at  the  cell  center,  for  example  j ,  E^ ^ ,  T^^ ,  e^.. , 
and  so  on.  The  horizontal  velocity  component  is  defined  at 
the  midpoint  of  vertical  cell  faces  and  the  vertical  component 
at  the  midpoint  of  horizontal  faces  (thus  we  have  u^+  ^ y  u^_ 

Vij+  hr  Vij-  J5) ' 

The  field  equations  (listed  in  section  Tl-f)  may  now  be 
written  in  finite-difference  form.  The  finite-difference  analogue 
of  the  mass  conservation  requirement  (11-53)  is: 

Preceding  page  blank 
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D.  . 
ID 


velocity  divergence  in  cell  ij 

_ 1 

r , 


I*?-  lri+%ui+%j-  ri-  Hui-  “ij1 


i  i 


+  Az j  ^Vij+  Jj"  Vi j  — *5 ^ 


(III-l) 


The  radial  and  vertical  mean  flow  momentum  equations  (11-54 
and  11-55)  become  respectively: 


(^i+  *  r^Ar.^  tri(uij>2-  ri+l(ui+lj)2] 


+  Az  j  Vi+Jsj-Js  ”Ui+  *s  Vi+,sj+  ^ 


^  [pii  '  pwi' 

_ L_  ,ri+lEi+lj  .  , 

ri+»s  ari+JiI  Ari+1  i+<*j  i+>ij 

* 

Ar^  ^i+Jsj  ”  ui-*ij^ 

*  * 

2ci+»sjui+^j  2 rfi+iili 

(1W2  Azj  Ari+Js 


(V  -  V  )- 

i+1  j  -*-*5  Vij+Js  Ari+Jj 


(vi+lj->s  '  Vij-J5)] 


+  Az  j  ^ei+Jsj+Js  (AZj+Jj  (ui+Jij+l_Ui+J5j  * 


(equation  continued  next  page) 
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r 


^  (vi+ij^  "  vij+%n 


i+Hi-h  ^ Az  j_js  (ui+%j  “  Ui+Jsj-1)_  Ari+}j  *vi+lj-Jj“vij-J*^  ] 


(HI-2) 


^3t^ij+*s  ri^ri  Ui-Jsj+J5  Vi_Jsj+J5  ri-*-J5Ui+3s j +JsVi+ *3 j '•-J5 ^ 


+  AZj+Js  Cvij  “  vij+l  1  +Azj+j5  tPij”Fij+l1+  ?ij+J59 


+  Az j+Jj  CAZj+1  (vij+|  ”vi j+Js^  ”  AzT1  (vij+Js”vij-Js)  1 


+  2  rri+^  E  (u  _u  ) 

+  r.Ar.  [  Az j+Jj  (Ui+>sj+l  Ui+J5j) 


ri-H  ei-hi+>i 


Azj^ 


(ui-Jsj+l  "Ui-Jsj)j 


+  riAri  tri+Jsei+Jsj+Js(Ari+Js  ( vi+l j+fs“Vi j-»*35) 


AZj+}j  ^ui+Jsj+l_ui+Jsj 


-  ri-*fi-y+%(Zr^  (vij^  ‘vi-lj^) 


AZj+Jj  (ui+»sj  +  l  "  Ui+J5j))j 


(HI-3) 


The  interpolated  quantities  appearing  in  the  above  equations 
are  defined  as  follows: 
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_  1  I  Az-i  +  i  Az-i-i  1  Azi 

Uij  *  Uij‘  2  5  )1  +  U,,x,I  ITT3—  ] 


‘j+<s  4rJ-* 


Lj+1  F  ^j+Jj 


+  u 


1 


ij-1  [F  AZj_^ 
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_ 1  ri+ri-J,  Ari-1  ri+k+r*  Ari 

Vii  =  F  vii  I  <•  r  *>  U+*  Tr^)  +  (~y- *>  (2-»S  T=i— ) 

l]  8  13  r±  Ari_ij  ri  Ari+jJ 


1  ri+ri-k  Ari-1 

+  F  vi-l-i 1  <-y— — >  <l-%  JT =i~i)  1 
8  1-I3 


i  r.  ,  ,+r.  ,  Ar. 

+  f  vi+ijt(— — ){  7  af: — )] 


'i+*s 


(III-5) 


where 


u 


=  .  (4^  +  1)])  +  u,  ,  .  + 


ij  ui+»jj'  r 


(III-6) 


Vij  *  2  ^ Vi j  +*s+vi j — *s  ^ 


(III-7) 


and 


_  Az.  Az. 

Ui+Jsj+Js  “  Ui+*ij  (1”*5  Az^J  +  ui+*s j+1  ^  Az^.  * 


j+*s 


*j+* 


(III-8) 
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(equation  continued 
next  page) 
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Ui+*5  j  +*5  “i+^j+Js 
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ij+Js 
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i+lj+^s  Ar 


i+*s 


(Ill-ll) 
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.  .*ri  *»i+l 


+  £i+lj  (4ri+)s4z j+)s 


+  CJ.x, ( 


aruiazi 
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*  *  Ari+1  *  Ari 

2£i+lij  -  +W^> 
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Azi+1  Azi 

j+Js  ^ij  (2Azj+Jj)+^ij+l  (iAZj+^5 


( III— 14 ) 


Furthermore,  forward  time  differencing  will  be  used  throughout; 
that  is,  for  example, 


N+l  N  N+Js  ,  3u.  N 

Ui+*j  "  Uij+Js  +  At  (Tt)i+>sj 


(III-15) 
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where  quantities  with  superscript  (N)  are  "old"  values  and 
with  (N+l)  "new"  values;  AtN+>*  is  the  time  interval  between. 


c.  Scalar  Transport 

The  scalar  transport  equations  are  handled  somewhat 
differently.  The  general  differential  form  for  scalar 
transport  is 


9Q 

at 


iL(r 

r  9r  1 


u  Q)  - 


a 

37 


(v°>  +  ?  If  <r*e  l?)+  k(YelF)+ 
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The  finite  difference  analogue  is  as  follows: 


(at*ij  (i-d.  .AtN+is)  IriAri  (Ci”,5Ui“,5]Qi",s3 


ri+Jsui+J5j®i+>sj  ^  +  A  ^vij-'s®ij-J5  vij+*s®i  j+*s^ 


+  r~Kr7  IAri+^yei+>s j  (Qi+lj_Qij) 


ri-H 
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ye 

Iz 
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and  we  will  again  update  values  using  forward  differences: 


(III-18) 


The  velocity  divergence  D^j  (see  III-l)  appears  in  the 
convection  terms  above  as  part  of  a  correction  factor.  In 


principle  D^j  should  always  be  identically  zero,  but  in 
the  calculational  scheme  this  is  not  always  true,  and 
small  velocity  divergences  may  appear,  as  will  be  discussed. 
This  correction  factor  tends  to  negate  convective  errors 
in  scalar  transport  arising  from  small  but  finite  velocity 
divergences. 

The  interpolated  concentrations  appearing  in  the 
convective  flux  terms  are  defined  as  follows: 

°U:  ■ 


* 

Qij+Js 


The  source  terms  for  the  various  scalar  quantities  are 
as  follows.  For  solutes,  there  are  no  sources,  so: 
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ns..  E  0 

ID 

For  turbulent  energy,  we  have  (see  11-56) 

2E,  ■ 
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ID  ID 
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and  for  temperature, 

2E . 

nT.  .  -  ?  ^  +  8'ij)  (III-22) 

ID  ij  J 

where  the  generalized  rate  of  mean  flow  shear  is  given  by: 
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d.  Turbulent  Mixing 

To  establish  the  eddy  viscosity  distribution  we  must 

first  fix  the  distribution  of  the  macroscale  A,  which  in 

turn  depends  on  certain  space  integrals  of  n2  (defined  above) 

and  the  quantity  (fifi')2as  discussed  in  section  II-c.  This 

latter  quantity  is  expressed  in  finite-difference  form  as 

follows:  (fl2  ft?  )! 

4(ftft')?. - i±±J - 

3  <AriV  4ri-l,>! 

(n2 . . , « —  n2 .  . ) 2 
+  i]-n  ip-r 

(4zj+^+Azj-is)  2 

The  "averaging  distance"  Ao  may  now  be  formed  as: 
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and  the  "characteristic  integrals"  I2  and  J2  ares 


all  cells 

Z  E  w.  .  ,  .  .  , 
i  •  j  •  1 1 1  i  J  t  D 

all  cells  ~  —  ™ 


Z  l  w 

i'  D' 


1/ 


2  2 
i  I  j  I  ) 


(III-26) 


Jij  = 


all  cells 

ij  j' _ 


all  cells 
Z  Z  w ,  .  ,  .  . , 
i  d 


m'2L, .,) 


(Ill— 27) 


where 
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and  4*(x)  is  defined  in  equation  (11-67);  it  is  illustrated 
in  Figure  3.  Now,  the  "scales  of  size"  may  be  formed. 

The  "microscale"  or  "dissipation  length"  is  just: 
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and  the  macroscale  is; 
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e.  Explicit  Solution  for  Pressure 

Examination  of  the  preceding  field  equations  shows 
that,  if  the  entire  set  of  velocities  and  Q's  is  known,  the 
changes  in  these  velocities  and  Q's  over  a  short  time 
interval  can  be  calculated.  The  only  missing  quantity  is 
the  pressure  distribution  (P^j).  In  incompressible  flow, 
the  pressure  distribution  may  be  found  from  the  mass  conserva¬ 
tion  constraint  (V*  u=  0).  We  first  write  an  expression 
for  the  time  rate  of  change  of  the  velocity  divergence  in 
cell  ij; 
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If  the  finite-difference  momentum  equations  (III-2  and  III-3) 
are  substituted  into  this  expression,  considerable  cancellation 
occurs,  and  the  result  is: 
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(Equation  continued  on  next  page) 
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The  pressure  equation  is  obtained  by  setting 

N+l 

(the  "new"  value)  equal  to  zero.  This  is  a 

"self -correcting"  feature  of  the  method;  even  if 

at  some  time  large  divergences  do  exist,  they  will 

be  quickly  eliminated.  That  is,  we  arrange  the 

rate  of  change  of  D. .  so  as  to  tend  to  "zero  out" 
N+l  1 J 

.  Therefore,  the  pressure  equation  is: 


ri±Js 


ri4ri 


“i-% 


,pii-pi-ii)I 


+  1 


(pij-pij+l>  - 


(PiJ-Pij-l»l 


+  Rj.  »  0  (III-35) 

where 


Note  that  contains  only  velocities,  eddy 
viscosities,  and  Sg's.  Therefore,  if  we  are  given 
a  set  of  these  variables,  equation  III-35  may  be 
solved  by  an  iterative  technique  for  the  pressure  field. 
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IV.  BOUNDARY  CONDITIONS 

a.  Fixed  Walls 

The  computing  mesh  is  bounded  within  an  upright  circular 
cylinder,  whose  walls  are  rigid  and  impermeable.  These  walls 
are  considered  to  be  "free-slip"  -  that  is,  along  a  wall,  the 
velocity  component  normal  to  the  wall  must  be  zero,  but  the 
tangential  component  is  not  constrained.  This  is  justified 
since  we  are  primarily  concerned  with  problems  at  high  Reynolds 
number,  and  it  is  to  be  expected  that  the  boundary  layer 
thickness  will  be  negligibly  small  compared  to  the  grid  resolution. 
This  "tank"  may  be  of  arbitrary  diameter  and  height,  so  that 
boundary  effects  may  be  minimized  by  moving  the  boundary  far 
away;  on  the  other  hand,  boundary  effects  may  be  of  interest 
(for  example,  the  effect  of  a  nearby  sea  bottom  upon  the  behavior 
of  an  underwater  explosion  bubble).  For  convenience,  we  will 
designate  the  boundaries  as  the  "wall",  "ceiling",  "floor", 
and  "axis".  The  axis,  of  course,  does  not  represent  a  true 
physical  boundary,  but  is  rather  an  axis  of  symmetry.  However, 
boundary  conditions  appropriate  for  a  physical  boundary  are  also 
appropriate  here. 

The  finite-difference  equations  when  applied  to  cell  ij  require 
quantities  located  in  neighboring  cells;  it  is  therefore  necessary 
to  add  one  layer  of  cells  outside  each  boundary  (see  Figure  4). 

The  boundary  conditions  are  then  applied  by  inserting  values  for 
the  field  variables  in  these  fictitious  cells  at  each  time  step 
such  that  the  boundary  conditions  are  satisfied.  The  field 
equations  provide  exactly  the  information  required  to  obtain 
these  relations.  For  simplicity  (and  without  loss  of  generality) 
we  use: 
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FIGURE  4:  LAYOUT  OF  THE  COMPUTING  MESH 


where  i  *  1  is  the  layer  "on  the  other  side"  of  the 
axis,  i  -  H  is  the  layer  outside  the  wall,  j  =  1 
is  the  layer  under  the  floor,  and  j  =  N  is  the  layer 
above  the  ceiling.  In  general,  the  boundary  conditions 
are  that  no  flux  of  transported  scalars  (whether 
convective  or  diffusive)  may  occur  through  a  boundary, 
and  that  the  velocity  constraint  mentioned  earlier 
applies.  These  will  be  satisfied  if  the  following 
obtains : 
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(Equation  continued) 
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(IV- 3) 
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where  Q  is  E,  T,  or  S. 


b.  Free  Surfaces 


The  boundary  conditions  at  a  free  surface  may  be 
summarized  as  follows: 

1)  There  must  be  no  stress  tangential  to  the  surface. 

2)  Pressure  must  be  continuous  across  the  surface. 

3)  There  must  be  no  Q  transport  through  the  surface. 
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The  problem  is,  where  are  these  conditions  to  be 
applied,  or  rather,  how  are  free  surfaces  to  be  located 
within  the  mesh?  In  order  to  identify  which  portion 
of  the  mesh  is  filled  with  fluid  and  which  portion  is 
empty,  we  insert  a  number  of  massless  "marker  particles" 
which  move  with  the  fluid  and  which  participate  in 
the  calculation  only  to  the  extent  that  they  allow  the 
computer  to  "flag"  cells  as  full  or  empty.  In  particular, 
at  the  end  of  a  time  step,  each  particle  is  moved  with  a 
velocity  which  is  interpolated  between  adjacent  principal 
velocity  points  in  the  mesh.  That  is,  designating  the 
coordinate  of  particle  k  as  r^,  z.  : 

For  r±  <  fk  <ri+1  and  2k<  *j+!i  . 


vk  =  vij->i 
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Then,  we  replace: 


+  v1 


+  V,-1 


A  tN+  ** 
itN+,i 


(IV-7) 


(IV-8) 


Two  sorts  of  void,  the  air  and  the  bubble,  are 
allowed  for  in  the  scheme,  and  therefore  we  have  two  free 
surfaces.  Cell  status  is  assigned  by  the  flagging  of  each 
cell  in  one  of  five  ways: 

AIR:  The  cell  is  empty  (contains  no  particles)  and 

is  in  the  region  designated  as  air. 

BUB:  The  cell  contains  no  particles  and  is  within 

the  bubble. 

AIRSUR:  The  cell  contains  particles,  but  is  directly 
adjacent  to  at  least  one  AIR  cell. 

BUBSUR:  The  cell  contains  particles,  but  is  directly 
adjacent  to  at  least  one  BUB  cell. 

FULL:  The  cell  contains  particles,  and  all  adjacent 
cells  also  contain  particles,  i.e.  are  either 
FULL,  AIRSUR,  or  BUBSUR. 

Note:  "Adjacent"  in  the  above  context  means  "sharing 

a  side  with"  -  diagonal  relationships  are  not 
considered. 

Thus,  at  the  beginning  of  a  problem,  each  cell  is  flagged 
in  one  of  the  above  ways.  At  the  end  of  each  time  cycle, 
the  particles  are  moved,  and  the  cells  reflagged  if  necessary. 

A  typical  arrangement  of  particles  and  cell  flags  is  shown 
in  Figure  5.  Only  certain  transitions  are  permitted,  however  - 
for  example,  a  cell  previously  flagged  as  FULL  may  not,  in  one 
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cycle,  change  to  an  AIR  cell.  This  rule  is  to  prevent  the 
spontaneous  opening  of  voids  deep  within  the  fluid  due  to 
the  happenstance  that  the  cell  should  momentarily  contain 
no  particles.  To  become  an  AIR  cell,  a  FULL  cell  must  first 
pass  through  AIRSUR  status.  These  transition  rules  are 
summarized  in  Figure  6. 

One  other  pathological  condition  can  occur,  however. 

A  cell  may  attempt  to  become  an  AIRSUR  and  a  BUBSUR  cell 
simultaneously  -  that  is,  it  contains  particles,  but  is 
faced  on  one  side  with  an  AIR  cell  and  on  the  other  by  a 
BUB  cell.  What  is  done  when,  and  if,  this  occurs  is  to 
assume  that  the  bubble  has  "leaked";  the  layer  of  fluid 
between  the  air  and  the  bubble  has  become  less  than  one 
cell  thick.  Computationally,  when  this  condition  is 
detected,  all  BUB  and  BUBSUR  cells  are  subsequently  re¬ 
designated  AIR  and  AIRSUR  cells  respectively.  It  should  be 
pointed  out,  in  passing,  that  there  is  no  requirement  that 
there  be  a  bvbble  involved  in  the  problem  at  all  -  the  MACYL 
scheme  can  handle  single-surface  problems  as  well  as  two- 
surface  ones. 

The  application  of  the  boundary  conditions  is  fairly 
straightforward,  once  the  flags  are  known  on  each  cell. 

For  pressure,  we  merely  say  that  the  pressure  in  an  AIR 
or  AIRSUR  cell  is  simply  PA  (an  input  constant)  and  that 
in  a  BUB  or  BUBSUR  cell  is  PD  (an  arbitrary  function  of 
bubble  volume  and/or  time) .  For  velocities  in  surface 
cells  of  either  type,  if  one  side  is  "open",  the  requirement 
that  =  0  supplies  the  missing  value.  If  more  than  one 
side  is  open,  we  require  that 

?  h  (ru)  *  0 
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FIGURE  5:  CELL  FLAGS 
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CELL  FLAG  TRANSITION  RULES 


V.  STABILITY  AND  ACCURACY 


As  is  true  of  finite-difference  schemes  in  general, 
the  MACYL  code  must  contend  with  potential  numerical 
instability.  These  instabilities  are  intrinsic  in  the  finite 
resolution  of  the  space-time  grid,  and  in  general  may  be 
found  by  a  comparison  of  the  finite-difference  equations 
with  the  original  differential  equations  using,  for  example, 
Taylor  expansions.  The  first  (and  most  obvious)  restriction 
is  the  incompressible  analogue  of  the  Courant  condition, 
which  arises  from  the  convection  terms  in  the  equations: 


At  <<  — —  everywhere  (V-l) 

|u| 

where  AX  is  a  space  interval,  and  u  is  a  velocity,  (either  a 
material  velocity  or  a  gravity  wave  phase  velocity) .  Thus  we 
obtain  two  requirements: 


where 

UMAX  =  raaxiraura  material  velocity 

W  =  longest  gravity  wavelength  present 

(usually  the  mesh  diameter) 

H  =  maximum  fluid  depth 
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Physically,  this  means  that  no  disturbance  is  allowed 
to  traverse  more  than  one  cell  in  a  single  time  step. 

Numerical  experiments  suggest  that  adequate  accuracy  can 
be  maintained  if  we  require  that: 

At  <  0.4  — -  ( V—  4 ) 

I  u  I 

1  'max 

Another  constraint  on  the  time  step  relates  to  rates 
of  diffusion.  To  prevent  instability,  we  must  require  that: 

At  <<  — -■  -  everywhere  (V-5) 

e* 

Actually,  this  requirement  will  only  be  important  when 
diffusive  effects  dominate  convective  effects:  for  the 
class  of  problems  of  interest  here,  (V-l)  will  generally 
override  (V-5) .  To  maximize  efficiency,  at  the  end  of 
each  time  cycle,  we  calculate  the  maximum  time  step  consistent 
with  numerical  stability  and  use  it  for  the  next  interval. 

Still  another  requirement  relates  to  the  value  of  the 
eddy  viscosity.  This  requirement  arises  from  high-order 
errors  in  the  convection  terms  of  the  finite-difference 
representation  of  the  momentum  equations.  These  high- 
order  errors  are  of  a  diffusive  character,  and  the  "artificial 
diffusivity"  thus  created  may  be  of  either  sign.  We  require 
that: 

e  >  |c(AX)2Q|  (V-6 ) 

where  c«0.7.  If  the  eddy  viscosity  at  a  particular  space- 
time  point  is  not  large  enough  to  satisfy  (V-6) ,  we  "boost" 
the  eddy  viscosity  so  as  to  just  satisfy  the  requirement. 

It  should  be  noted,  however,  that  this  "boosting"  occurs 
only  in  the  "vorticity-dif fusion"  terms  of  the  momentum 
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equation.  That  is,  the  turbulent-plus-viscous  terms 
in  the  momentum  equations  are  of  the  form: 
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-f  1  3  ,  *  3u. 

21  F  5F  (r£  5F>  -  r 


* 

eu 


±  (e*  iZ)l 

3  Z  1  3  Z '  J 


+ 


* 

e  ' 


,3u  3v.  , 

3 z  "  3r  J 


z-direction: 


2[ 


13  .  *  3u, 

r  IF  (r£  IF> 


* 

(e 


3 

3r 


[r  e*' 


,3u  3vx , 
3z  "  3rn 


and  the  boosted  viscosity  only  appears  in  the  finite- 

difference  analogues  of  the  equations  where  indicated  by 

* 

e  ' .  Violations  of  this  requirement  tend  to  generate 
spurious  wave-like  disturbances  or  vortices  of  wavelength 
one  computational  cell. 

One  test  of  the  accuracy  of  the  method  is  whether  or 

not  transported  quantities  (momentum,  salinity,  etc.)  are 

conserved  rigorously.  It  may  be  shown  directly  from  the 

finite-difference  equations  that  this  is  so,  and  therefore 

-15 

the  error  is  of  order  computer  round-off  error  (about  10 
Another  test  of  any  incompressible  scheme  is  whether  mass, 
or  fluid  volume,  is  conserved.  If  were  always  exactly 
zero,  this  would  be  true  by  definition.  As  has  been  shown 
however,  it  is  not  generally  zero,  but  only  of  very  small 
magnitude  (D^j  may,  of  course,  be  made  arbitrarily  small 
by  "tightening"  the  convergence  criterion  in  the  pressure 
iteration  procedure)  .  A.nother  way  of  checking  on  mass 
conservation  is  to  keep  track  of  the  total  fluid  volume  as 
the  calculation  proceeds.  There  is  a  slight  ambiguity 
here,  however.  The  problem  is,  for  a  surface  cell,  how 
much  of  the  cell  is  considered  to  be  full  of  fluid? 
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Somewhat  arbitrarily,  the  following  rules  were  chosen  to 
evaluate  the  best  estimate  of  the  fluid  volume  in  a 
surface  cell: 

1)  If  the  cell  is  open  on  one  side  only,  or  on  two 
opposite  sides,  it  is  considered  to  be  half-full. 

2)  Otherwise,  it  is  considered  one-quarter  full. 
Naturally,  total  volumes  calculated  in  this  way  fluctuate 
slightly  from  cycle  to  cycle  due  to  the  surface  ambiguity. 

The  amplitude  of  the  fluctuation  decreases,  of  course,  as 
resolution  is  increased.  It  has  been  shown,  however,  that 
even  for  problems  carried  out  through  thousands  of  time  cycles, 
the  mean  of  the  fluctuating  values  does  not  shift,  but 
remains  the  same  throughout  the  calculation.  Therefore, 

the  small  D^'s  are  in  some  sense  randomly  distributed,  and 
their  effects  cancel. 
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VI.  THE  COMPUTER  PROGRAM 


The  necessary  tools  have  now  been  developed  to  assemble  a 
procedure  for  solution.  If,  at  time  t,  we  know  the  entire 
state  of  the  system,  we  may  update  the  field  variables  to  their 
values  at  t+  At  through  the  following  procedure: 

Step  1  -  Calculate  and  store  values  of  R^ ^  (defined  by  equations 
III-34  HI-36)  for  all  FULL  cells. 

Step  2  -  Iterate  a  solution  for  the  entire  pressure  field  in 

FULL  cells  using  a  Gauss-Seidel  iteration  procedure  in 
connection  with  equation  III-35  until  the  entire  P^ 
field  converges.  The  criterion  for  convergence  is: 


Step  3  - 
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for  all  ij 

where  superscripts  k-1,  k  denote  successive  passes 
through  the  iteration  loop,  and  H  is  the  maximum  fluid 
depth  expected  in  the  problem. 

Calculate  new  velocities  based  on  the  new  set  of 
pressures,  the  momentum  equations  (III-2  and  -3),  and 
the  boundary  conditions  discussed  in  section  IV. 
Calculate  and  store  values  for  the  quantities 
and  for  all  FULL,  AIRSUR  and  BUBSUR  cells 

(IH-23,-24)  . 

Compute  the  "averaging  distance"  Ao  per  equation  III-25. 
Calculate  the  "macroscale"  (A^j)  and  "microscale" 

(A^j)  distributions  for  all  FULL,  AIRSUR,  and  BUBSUR 
cells  as  outlined  in  equations  III-26  through  III-30. 
evaluate  the  eddy  viscosity  distribution  (e^j)  in  all 
FULL,  AIRSUR,  and  BUBSUR  cells  (III-31). 
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Step  8 


Step  9 
Htep  10 
Step  11 


Step  12 


Utep  13 
Step  14 


Update  time  (tN+^  *  tN  +  A  tN+J*)  ;  calculate  a  new 
time  interval  based  on  stability  requirements  V-4 
and  V-5;  output  the  state  of  the  system. 

Update  E^j  (the  turbulent  energy  distribution) 
using  equations  III-17  and  III-21. 

Update  the  temperature  distribution  (the  T\j's) 
using  equations  III-17  and  III-22. 

Update  solute  concentrations  using  III-17  and  III-20. 
The  program  allows  for  two  solute  fields;  one  is 
usually  salinity,  and  the  other  may  be  used  to 
follow  a  dissolved  contaminant. 

Re-evaluate  ^ j  for  all  FULL,  AIRSUR  and  BUBSUR 
cells  using  the  new  T^j's  and  S^^'s,  and  the 
appropriate  "equation  of  state"  for  the  fluid 
(such  as  equation  11-52  for  seawater) . 

Move  the  marker  particles  according  to  IV-5,  -6, 

-7,  and  -8. 

Reflag  cells  as  required,  as  discussed  in  section 
IV-b;  re-evaluate  P  (the  bubble  pressure)  and  return 
to  step  1. 
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VII.  CONCLUSIONS 


The  MACYL6  computer  program  is  currently  operating  on 
the  CDC  6600  computing  system  at  the  Lawrence  Radiation 
Laboratory  in  Berkeley,  California.  As  has  been  discussed, 
earlier  versions  of  the  scheme  have  been  checked  against 
both  analytic  solutions  and  experimental  data,  with  favorable 
results;  in  general,  the  agreement  with  data  is  within  error 
of  measurement.  The  program  will  shortly  be  put  to  work  in  an 
investigation  of  the  explosion  debris  transport  from  very 
deep  underwater  nuclear  explosions.  The  potential  of  the 
scheme  for  application  to  other  problems  associated  with  under¬ 
water  explosions,  as  well  as  oceanographic  and  meteorological 
problems  in  general,  appears  to  be  great. 
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APPENDIX  -  LIST  OF  SYMBOLS 


C  molecular  diffusion  coefficient 
velocity  divergence  in  cell  ij 
E  Turbulent  kinetic  energy  per  unit  mass 


e 

-*■ 

g 

g 

H 


i,  j,k 


i»  j 


I 

J 

M 

N 

N 


2 

2 


,i:+i 


2.71828.  .  . 

acceleration  of  gravity  (vector) 

vertical  component  of  acceleration  of  gravity 

maximum  fluid  depth  expected  in  problem 

Part  II:  Cartesian  tensor  indices  (subscripts) 

Part  III:  et  seq. :  space  grid  indices  in  r,  z  directions 
respectively  (subscripts) 

first  characteristic  integral 

second  characteristic  integral 

maximum  value  of  i  index  (denotes  "wall")  (subscript) 
maximum  value  of  j  index  (denotes  "ceiling")  (subscript) 

integer  index  denoting  time  step  (superscript) 
mean  total  pressure  =  $  +  §E 

air  pressure 
bubble  pressure 

instantaneous  generalized  scalar  concentration 

mean  generalized  scalar  concentration 

turbulent  generalized  scalar  concentration  fluctuation 

'ij  +  Dij/4t 

radial  coordinate 

radial  coordinate  of  center  of  cell  ij 
radial  coordinate  of  outer  boundary  of  cell  ij 


?k  radial  coordinate  of  particle  k 


S'  instantaneous  salinity  (or  solute  concentration) 
S  mean  salinity  (or  solute  concentration) 
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e+  c+  c+  «+ 


T' 

T 


i 


u 


W 

w 


AX 

z 


3 

Zj+H 

~zk 


instantaneous  temperature 

mean  temperature 

time 

instantaneous  velocity  (vector) 
mean  velocity  (vector) 

turbulent  velocity  fluctuation  (vector) 

radial  velocity 

radial  velocity  of  particle  k 

vertical  velocity 

vertical  velocity  of  particle  k 

longest  gravity  wavelength  present 

weighting  function 

generalized  space  interval 

vertical  coordinate 

vertical  coordinate  of  center  of  cell  ij 
vertical  coordinate  of  upper  boundary  of  cell 
vertical  coordinate  of  particle  k 


ij 


or 

3 

Y. 


3 

A 

V 


e 

e* 


K 


turbulence  model  functions 

mean  flow  strain  rate  tensor 
Kroeneker  delta  function 
denotes  partial  differentiation 
denotes  a  finite  difference  or  interval 
vector  nabla  operator 
eddy  viscosity 
£  +  U 

molecular  thermal  diffusivity 
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! 


A0 ,Ai , .. successive  approximations  to  A 
A  turbulent  macroscale 

X  turbulent  microscale 

u  kinematic  viscosity 

5  density  variation  factor  ■  1  + 

IIq  source  term  for  scalar  quantity  Q 

tt  3.14159265... 

p  fluid  density 

o  specific  heat 

<t>'  instantaneous  pressure 
<!>  mean  pressure 

<)>  turbulent  pressure  fluctuation 

Y  geometric  function  for  cylindrical  coordinates 

ft  generalized  mean  flow  strain  rate 

ft'  generalized  mean  flow  strain  rate  gradient 

i 
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ERRATA 


IRA-TR-1-70,  "The  MACYL6  Hydrodynamic  Code:  A  Numerical 
Method  for  Calculating  Incompressible  Axisymmetric 
Time-Dependent  Free-Surface  Fluid  Flows  at  High  Reynolds 
Number,"  by  John  W.  Pritchett,  15  May  1970. 


p.  9:  sentence  after  Eq.  (11-27)  should  read: 

"Difficulty  will  of  course  be  experienced  in  calculating 
the  A  distribution,..." 

p.  12:  Eq.  (11-43)  should  read: 


